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ABSTRACT 


A  theoretical  study  of  the  restricted  primitive  model  of  the  electrical  double  layer 
using  a  free  energy  density  functional  theory  is  presented.  The  ion-ion  hard-core  repulsive 
contribution  to  the  free  energy  is  incorporated  by  a  nonlocal  excess  hard-sphere  free  energy 
density  functional.  The  electrostatic  part  of  the  ion-ion  direct  correlation  function  of  the 
inhomogeneous  electrolyte  in  the  interfacial  region  is  approximated  by  that  of  the  homo¬ 
geneous  bulk  electrolyte.  The  generalized  van  der  Waals  model  and  the  Tarazona  model 
are  used  to  construct  the  respective  excess  hard-sphere  free  energy  density  functionals 
and  compared  with  each  other.  Each  model  requires  as  input  a  hard-sphere  equation  of 
state.  The  Camahan-Starling  equation  of  state  is  chosen.  The  theory  predicts  correctly 
the  layering  effect  of  the  counterions  and  the  charge  inversion  phenomenon.  The  results 
for  1:1  and  2:2  electrolytes  agree  well  with  the  Monte  Carlo  simulations.  The  Tarazona 
model  gives  results  in  closer  agreement  with  Monte  Carlo  data  than  does  the  generalized 
van  der  Waals  model.  The  diffuse  layer  potentials  predicted  by  the  Tarazona  model  are 
not  as  accurate  as  those  by  the  generalized  hard-rod  model  which  has  been  applied  to  the 
same  problem  in  a  previous  study. 


I.  INTRODUCTION 


The  electrical  double  layer  is  the  inhomogeneous  region  of  an  electrolyte  near  a 
charged  surface.  Its  most  distinguishing  feature  is  the  charge  separation  in  the  electrolyte 
in  response  to  the  electric  field  produced  by  the  surface  charges.  Understanding  the  elec¬ 
trical  double  layer  is  fundamental  to  a  variety  of  phenomena  in  electrochemical,  colloid, 
and  biological  sciences. 

Because  of  many  complications  present  in  physical  systems  containing  the  electrical 
double  layers,  theoretical  studies  have  concentrated  largely  on  some  simplified  models.  The 
advantage  in  using  a  simplified  model  is  that  it  allows  an  unambiguous  test  of  a  theory 
against  the  computer  simulations  for  the  same  model.  Among  these  models,  the  most 
widely  studied  one  is  the  restricted  primitive  model  (RPM),  in  which  the  electrolyte  is 
represented  by  a  fluid  of  charged  hard  spheres  with  diameter  d  and  point  charges  qa  at 
their  centers,  the  solvent  by  a  continuum  with  an  isotropic  dielectric  constant  e,  and  the 
electrode  by  an  infinite  planar  hard  wall  with  a  uniform  surface  charge  density  a.  The 
electrode  surface  is  also  assumed  to  be  impenetrable.  The  RPM  of  the  electrical  double 
layer  is  completely  defined  once  the  system  parameters  e,  a,  d ,  {<?<*},  the  ion  number 
densities  in  the  bulk  electrolyte  {n°  },  and  the  temperature  T  are  specified. 

Studies  of  the  RPM  of  the  electrical  double  layer  began  with  the  pioneering  works  of 
Gouy  [1]  and  Chapman  [2]  based  on  the  Poisson-Boltzmann  equation,  in  which  the  ions 
were  treated  as  point  charges  so  that  the  hard-core  repulsions  between  the  ions  and  between 
the  ion  and  the  charged  wadi  were  not  included.  A  modification  was  made  later  by  Stem 
[3]  to  incorporate  the  ionic  size  effects  in  such  a  way  that  sill  ions  were  allowed  to  approach 
only  up  to  a  closest  distance  from  the  charged  surface,  whereas  the  ion-ion  interactions 
were  still  treated  by  the  point  charge  model.  This  was  called  the  modified  Gouy-Chapman 
(MGC)  theory  [1-3].  The  MGC  theory  is  quantitatively  accurate  for  weakly  coupled  (say, 
small  a  and  qa)  systems  at  low  densities,  in  which  case  the  effects  of  finite  ion  sizes  are  not 
important,  but  fails  to  represent,  even  qualitatively,  highly  coupled  systems  or  high  density 
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ones.  Qualitatively,  there  are  two  fundamentally  important  phenomena  present  in  the  ion 
density  profiles  of  the  RPM  of  the  electrical  double  layer  as  revealed  by  the  Monte  Carlo 
(MC)  simulations  of  Torrie  and  Valleau  [4]:  the  layering  effect  of  the  counterions,  and  the 
charge  inversion.  Firstly,  for  1:1  electrolytes  near  a  highly  charged  surface  the  counterion 
density  profile  is  no  longer  monotone  but  shows  a  second  layer  next  to  the  inner  monolayer. 
Secondly,  for  1:1  and  2:2  electrolytes  at  high  concentrations  and  high  surface  charges  there 
is  charge  inversion,  i.e.,  the  coion  density  exceeds  the  counterion  density  within  certain 
parts  of  the  interfacial  region.  The  MGC  theory,  however,  always  predicts  monotonic 
profiles  of  the  ion  density  [5].  Accordingly,  various  efforts  have  been  directed  toward 
interpreting  of  the  layering  and  the  charge  inversion  effects  through  better  treatments  of 
both  the  short-ranged  hard-sphere  interaction  and  the  long-ranged  Coulombic  interaction 
in  the  electrical  double  layer  [5,6].  These  efforts  include  works  based  on  the  modified 
Poisson-Boltzmann  (MPB),  the  Yvon- Born-Green  (YBG),  and  the  Omstein-Zernike  (OZ) 
equations.  There  are  also  works  based  on  free  energy  density  functional  theories.  We 
sketch  in  Table  I  these  different  theories  of  the  RPM  of  the  electrical  double  layer.  For 
detailed  references,  see  the  recent  reviews  [5,6]  of  this  subject  and  a  recent  paper  by  some 
of  the  present  authors  [19], 

The  MPB  theory  is  based  on  the  Kirkwood  hierarchy  and  requires  a  closure  which 
connects  the  density  distribution  functions  and  the  mean  electrostatic  potential  function. 
Computations  based  on  this  theory  have  proved  efficient  [5].  The  best  and  most  recent 
version  of  the  MPB  theory  is  the  MPB5  theory  of  Outhwaite  and  Bhuiyan  [7],  which  is 
successful  for  1:1  electrolytes  but  not  for  2:1  and  2:2  electrolytes  at  similar  densities  and 
temperatures.  The  YBG  theory  has  a  conspicuous  feature:  it  satisfies  the  contact  theorem 
exactly.  For  closure,  this  theory  requires  approximations  to  the  pair  distribution  function. 
A  promising  nonlocal  closure  has  been  proposed  recently  [8].  The  results  of  the  YBG 
theory  are  good  at  high  concentrations  but  poor  for  dilute  solutions;  this  is  believed  to 
result  from  numerical  difficulties  when  the  density  profiles  become  very  long-ranged  in  the 
low  concentration  limit  [8].  The  OZ  equation  together  with  the  mean  spherical  approxi¬ 
mation  (MSA)  or  the  hypernetted-chain  (HNC)  approximation  underlies  the  MSA/MSA. 


4 


the  HNC/MSA,  and  the  HNC/HNC  theories,  of  which  the  HNC/MSA  is  the  best.  In  the 
various  versions  of  this  theory,  the  ion-ion  direct  correlation  function  of  the  inhomogeneous 
electrolyte  in  the  interface  is  approximated  by  that  of  the  homogeneous  bulk  electrolyte. 
Numerical  predictions  from  the  HNC  theories  require  large  computational  efforts,  but  these 
theories  have  the  advantage  that  they  can  be  generalized  easily. 

Two  recent  works  on  the  extensions  of  the  HNC  theories  of  the  RPM  of  the  electrical 
double  layer  warrant  mention.  One  is  the  local  HNC  (HNCL)  and  the  local  MSA  (MSAL) 
approximations  suggested  by  Forstmann  et  al.  [9],  who  derived  an  integral  equation  for  the 
density  distribution  functions  from  a  density  functional  theory,  in  which  the  nonuniform 
Omstein-Zemike  direct  correlation  function  should  be  approximated.  The  essence  of  these 
approximations  is  that  the  exact  ion-ion  direct  correlation  function  of  the  inhomogeneous 
electrolytes  is  replaced  locally  by  either  the  HNC  or  the  MSA  solutions  to  the  direct 
correlation  functions  of  some  homogeneous  non-neutral  electrolytes.  These  non- neutral 
electrolytes  are  chosen  in  such  a  way  that  they  have  the  ion  compositions  corresponding 
respectively  to  those  found  locally  in  different  regions  of  the  double  layer  [9].  These  theories 
give  good  results,  but  an  unsatisfactory  feature  of  these  approaches  is  that  they  have 
one  undetermined  parameter  which  must  be  fitted  to  the  MC  simulations.  The  other 
recent  extension  of  the  HNC  theories  is  that  of  Plischke  and  Henderson  [10],  who  solved 
the  inhomogeneous  OZ  equation  in  the  HNC  and  the  MSA  closures  simultaneously  with 
the  Lovett-Mou-Buff-Wertheim  equation  for  the  ionic  density  profiles  (OZ/LMBW).  The 
OZ/LMBW  theory  gives  the  most  accurate  predictions  currently  available.  However,  its 
numerical  implementation  involves  rather  large  computing  time  and  considerable  storage. 
In  summary,  all  the  theories  mentioned  above  succeed  quite  well  for  weakly  coupled  systems 
at  low  densities,  whereas  only  the  YBG,  the  HNCL  and  MSAL,  and  the  OZ/LMBW 
theories  predict  correctly  the  qualitative  behavior  of  highly  coupled  systems  at  moderate 
or  high  densities.  Quantitative  agreement  with  the  MC  results  is  accomplished  without 
adjustable  parameters  only  by  the  OZ/LMBW  theory. 

In  contrast  with  the  approaches  mentioned  above,  there  are  far  fewer  investigations 
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based  on  the  density  functional  theories.  The  essence  of  the  density  functional  theory  is  to 

construct  the  Helmholtz  free  energy  as  a  functional  of  the  particle  density  and  minimize 

the  free  energy  to  obtain  an  expression  for  the  equilibrium  density  distribution.  In  the 
* 

currently  tractable  versions  of  this  approach,  the  contribution  of  particle-particle  pair 
interactions  to  the  free  energy  are  separated  into  hard-sphere  and  electrostatic  parts.  The 
hard-core  repulsion  between  the  ions  is  represented  by  a  nonlocal  generic  free  energy'  density- 
functional  for  inhomogeneous  hard-sphere  fluids,  which  was  first  proposed  by  Percus  [11] 
and  generalized  to  fluid  mixtures  by  Vanderlick  et  al.  [12].  Several  models  of  such  a 
nonlocal  free  energy  density  functional  have  evolved  in  recent  years.  These  include  the 
generalized  van  der  Waals  (GVDW)  model  [13],  the  Tarazona  (TRZ)  model  [14],  and  the 
generalized  hard-rod  model  (GHRM)  [15,16].  Each  of  these  models  requires  as  input  a 
hard-sphere  equation  of  state.  A  comparison  of  these  models  with  applications  to  hard 
sphere  fluids  confined  between  hard  walls  has  been  reported  [17].  The  first  application  of 
the  density  functional  theory  to  the  RPM  of  the  electrical  double  layer  was  the  work  by 
Boyle  et  al.  [18],  in  which  the  GVDW  model  together  with  the  Clausius  equation  of  state 
was  used  and  the  electrostatic  interaction  between  the  ions  was  treated  in  the  mean  field 
approximation.  The  results  obtained  compare  somewhat  poorly  with  the  MC  simulations. 
This  can  be  attributed  in  part  to  two  facts:  the  Clausius  equation  of  state  does  not  account 
well  for  the  behavior  of  a  hard  sphere  fluid;  and  the  mean  field  approximation  is  too  crude 
to  represent  accurately  the  electrostatic  correlations.  Recently,  Mier-y-Teran  et  al.  [19] 
have  applied  the  GHRM  together  with  the  Caroahan-Starling  equation  of  state  to  the  same 
problem.  An  improved  treatment  of  the  electrostatic  interaction  was  used,  in  which  the 
electrostatic  part  of  the  ion-ion  direct  correlation  function  of  the  inhomogeneous  electrolyte 
in  the  interfacial  region  was  approximated  by  that  of  the  homogeneous  bulk  electrolyte. 
This  electrostatic  correlation  consists  of  a  long-ranged  mean  field  part  and  a  short-ranged 
residual  part.  Also,  the  volume  averaging  technique  of  Ref.  18  as  it  applied  to  regions  very 
close  to  the  charged  surface,  which  is  now  believed  inappropriate,  has  been  eliminated. 
The  GHRM  has  am  accuracy  approaching  in  some  cases  that  of  the  OZ/LMBW  theory, 
which  is  the  best  available,  but  requires  only  a  small  fraction  of  the  computing  time  used 
to  solve  the  OZ/LMBW  theory.  This  is  in  general  an  attractive  feature  of  the  density 


functional  theory,  especially  in  the  light  of  the  eventual  goal  of  treating  the  solvent  as  a 
molecular  fluid  instead  of  a  dielectric  continuum. 

It  has  been  shown  that  the  TRZ  model  is.  overall,  a  better  approximation  for  the 
inhomogeneous  hard-sphere  fluids  than  the  GHRM  [17],  and  that  the  GHRM  predicts 
qualitatively  correct  results  for  the  RPM  of  the  electrical  double  layer  although  there 
are  small  quantitative  differences  from  the  MC  data  [19].  A  next  logical  step  in  the 
development  of  density  functional  theory  is  to  examine  the  accuracy  of  the  TRZ  model  in 
applications  to  the  RPM  of  the  electrical  double  layer.  The  purpose  of  this  paper  is  to 
carry  out  such  an  examination. 

In  section  II,  we  use  the  nonlocal  free  energy  density  functional  theory  to  derive  an 
integral  equation  for  the  ion  density  distribution  functions  of  the  RPM  of  the  electrical 
double  layer.  We  extend  the  TRZ  model  to  mixtures  of  hard  spheres  of  the  same  size  and 
then  apply  it  together  with  the  Carnahan-Starling  equation  of  state  to  construct  the  excess 
hard-sphere  free  energy  functional.  Because  of  its  simplicity,  the  GVDW  model,  which  is 
in  fact  the  lowest  order  approximation  to  the  TRZ  model,  is  also  examined.  The  approach 
we  use  is  the  same  as  that  in  Ref.  19.  The  results  for  the  ion  density  profiles  and  the 
mean  electrostatic  potential  profiles  are  reported  in  Sec.  Ill,  and  compared  with  the  MC 
simulations  [4]  and  the  results  of  other  theories,  especially  those  of  Ref.  19.  In  Sec.  IV.  we 
conclude  with  a  brief  summary  and  suggest  improvements  on  the  present  approach.  We 
also  mention  the  extensions  of  the  density  functional  theory  to  more  realistic  models  of 
the  electrical  double  layer. 


II.  THEORY 


A  general  starting  point  for  the  treatment  of  the  inhomogeneous  fluids  is  the  density 
functional  theory  in  which  the  Helmholtz  free  energy  F  is  expressed  as  a  functional  of  the 
particle  density  distribution  n(r)  [20,21].  For  a  system  of  charged  particles  of  different 
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kinds  a  (a  =  1,2, . . .  ,m)  at  fixed  temperature  T,  volume,  the  external  field  of  potentials 
ra(r),  and  the  chemical  potentials  /iQ,  the  equilibrium  particle  density  distribution  is  that 
which  minimizes  the  grand  potential  17  of  the  system  where 

i 

17({n})  =  F({n})  -  ]T  J  d3rnana{ r).  (2.1) 


The  equilibrium  particle  density  distribution  is  thus  a  solution  to  the  following  Euler- 
Lagrange  equation  resulting  from  the  minimization  of  17({n})  with  respect  to  na{ r): 


6 17( {n})  _  <$F({n}) 
Sna{ r)  ”  6na( r) 


(2.2) 


The  Helmholtz  free  energy  density  functional  of  the  system  mentioned  above  is  as¬ 
sumed  to  be' written  as 


F({n})  =  [  d3r  n«(  r)  [/n  (A*  na(r))  -  l] 

a  J 

+  5^  J d3ma( r)wo(r)  -  <^({n}), 

a  J 


(2.3) 


where  k  is  Boltzmann’s  constant  and  Aa  denotes  the  de  Broglie  wavelength  of  particles 
of  species  a.  The  first  term  on  the  right  side  of  Eq.  (2.3)  is  the  ideal  gas  contribution 
to  the  free  energy.  The  second  term  corresponds  to  the  external  field  contribution.  The 
third  term,  the  functional  —  <£({n}),  represents  the  contributions  from  the  ion-ion  pair  in¬ 
teractions.  The  functional  ^({n})  is  a  unique  functional  of  particle  density  in  a  specific 
system  and  its  derivatives  with  respect  to  density  define  a  hierarchy  of  direct  correla¬ 
tion  functions  [20,21].  Of  particular  importance,  the  nonuniform  Omstein-Zemike  direct 
correlation  function  cQg( r,  r7)  is  given  by 


Caf)( r,  r') 


1  <$2^({n}) 

kT  6na(r)6n^{r')' 


(2.4) 


If  we  have  some  means  of  evaluating  the  direct  correlation  functions  over  a  range  of  nonuni¬ 
form  densities,  we  cam  obtain  <£({n})  by  functional  integrations  of  the  direct  correlation 
functions  with  respect  to  density. 
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Although  <?>(  { n } )  does  not  depend  on  the  specific  choice  of  integration  path  [20],  it 
is  most  convenient  to  integrate  along  a  linear  density  path  from  some  initial  density  n^(r) 
to  the  final  density  riQ(r), 

na(r;  A)  =  n^(r)  +  A  fna(r)  -  n'a{ r)]  ,  (2.5) 


where  A  is  a  parameter  varying  between  0  and  1.  Such  an  integration  can  be  further 
simplified  if  the  initial  density  is  chosen  to  be  zero,  n^r)  =  0.  Because  the  particle-particle 
pair  interactions  become  negligible  at  sufficiently  low  density,  we  have  <s>(  { n’ } )  =  0  in  this 
case  [21].  This  yields 


0({n})  =  fcTV  /  f  d3rd3r'na(r)ng(r')  f  d\  fdA'c^r.r';  A'). 


>.6) 


A  reference  (initial)  state  different  from  that  used  in  Ref.  19  has  been  chosen  here  in 
the  expression  of  the  density  functional  but  the  equivalence  of  the  two  approaches  is 
guaranteed  by  the  uniqueness  of  <t>. 


The  free  energy  density  functional  then  becomes 
F({n})  =  kT  ?/  d3rna( r)  [/n  (Aana(r))  -  l] 

+  ?/  d3r  nQ(r)ua(r) 

J  d3rd3r'na(r)ng(r')  dX  d\' cag(r,  r';  A'). 


(2.7) 


It  bears  emphasizing  that  in  this  paper  we  employ  the  Helmholtz  free  energy  density 
functional  rather  than  the  grand  potential  functional  which  was  used  in  Ref.  19.  These 
two  functionals  are  simply  related  to  each  other  by  Eq.  (2.1).  No  approximations  have 
been  made  so  far  concerning  the  physical  system  we  are  going  to  investigate,  and  in  this 
sense  Eq.  (2.7)  is  exact. 


For  an  inhomogeneous  RPM  electrolyte  with  both  hard-core  repulsions  and  Coulom- 
bic  interactions,  the  nonuniform  Orastein-Zemike  direct  correlation  function  cog(r,  r')  is 
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not  in  general  known.  But  it  cam  be  decomposed  into  three  parts: 

(2-8) 

The  first  term  on  the  right  side  of  Eq.(2.8),  c^(r,  r'),  is  the  direct  correlation  function 
of  a  fluid  of  neutral  hard  spheres  due  solely  to  haxd-core  repulsion,  which  is  discussed  in 
more  detail  in  what  follows.  The  second  term  represents  the  ionic  correlation  due  solely 
to  the  Coulombic  interaction  between  charged  particles.  Here  we  have  made  an  additional 
assumption  about  the  RPM  of  the  electrical  double  layer:  the  electrode  surface  is  taken 
as  having  the  same  dielectric  constant  e  as  has  the  solvent.  In  this  case  there  are  no  im¬ 
age  effects  arising  from  the  polarizations  of  the  surface,  and  the  electrostatic  interaction 
between  ions  remains  in  the  form  of  central  force  as  expressed  by  Eq.(2.8).  An  additional 
contribution  in  the  same  equation,  Ac0(j(r,r'),  can  be  interpreted  in  some  sense  as  the  re¬ 
sults  of  the  cross- correlations  between  the  Coulombic  and  the  hard-sphere  interactions.  To 
proceed  further  we  must  make  some  kinds  of  approximations  for  c^(r,  r')  and  &co9(r.  r'). 
We  choose  here  to  approximate  Aca^(r,  r1)  of  the  inhomogeneous  electrolyte  in  the  inter¬ 
face  by  that  of  the  homogeneous  neutral  bulk  electrolyte,  AcQ/s(|r  —  r'|).  The  latter  is 
given  by  the  analytical  solution  to  the  OZ  equation  when  the  MSA  closure  is  used  [22]. 
The  explicit  expression  for  AcQ^(|r  —  r'|)  is  given  in  Ref.  19. 

On  the  other  hand,  instead  of  making  any  kind  of  explicit  approximation  to  the 
hard-sphere  correlation  function  c£^(r,  r'),  we  introduce  an  excess  free  energy  density 
functional  Fezcett({a})  into  our  formalism  to  incorporate  the  hard-core  repulsion  between 
the  ions.  F**ce**({n})  is  essentially  the  excess  free  energy  of  a  hard  sphere  fluid  over 
the  ideal  gas  contribution  accounting  for  the  finite  size  of  the  particles.  This  excess  part 
of  the  free  energy  is  extremly  important  in  determining  the  structure  of  Coulombic  and 
non-Coulombic  inhomogeneous  fluids  provided  there  is  short-ranged  hard-core  repulsion 
between  the  particles. 

On  the  basis  of  all  the  foregoing,  the  free  energy  density  functional  can  be  rewritten 


cag( r,r')  =  c^(r,r')  - 


Qaq$ 


e&r|r  —  r'| 


+  -^cQ^(r,  r'). 
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F(<n})  =  ir£ 

a 

'  +  f  d3rn°‘(r)vM 

a 

+  ^^3r'n0(f)nj(r')^ 7| 

a/J 

J  J  d3rd3r'na(r)nd(r')Aca3(\r -r'\) 

a$ 


j  d3rna( r)  [/n  (A’na(r))  -  l] 


(2.9) 


+  FeIcejJ({n}). 

Here  the  external  field  potential  va(r)  in  the  planar  electrical  double  layer  can  be  shown 
to  be  [5] 

va(r)  =  -2nqacrx/e  +  t?*'(x),  (2.10) 


where  x  =  |r  •  x|  is  the  ion’s  distance  to  the  wall, 

vks(x)  =  oo,  x  <  d/2, 

=  0,  x  >  d/2. 


(2-11) 


As  an  intuitive  generalization  of  the  exact  theory  of  the  inhomogeneous  hard-rod 
fluids,  the  excess  hard-sphere  free  energy  functional  Feice"({n})  suggested  by  Percus  is  in 
general  a  nonlocal  density  functional  and  can  be  expressed  in  terms  of  two  coarse-grained 
densities,  n£(r)  and  ri£(r),  which  are  themselves  at  each  point  density  functionals  of  the 
local  density  distribution  na( r): 

F“c’"({n})  =  E/ *'WrVU*T(r))'  (212) 

where  Fo(n)  is  the  excess  free  energy  per  particle  of  a  homogeneous  hard  sphere  fluid  of 
density  n,  and  its  expression  can  be  derived  from  the  hard-sphere  equation  of  state  [10,11]. 
The  physical  interpretation  of  Eq.(2.12)  is  that  the  excess  free  energy  contributed  from 
a  unit  volume  around  r  is  the  product  of  two  quantities:  the  average  number  density  of 
particles,  nv( r),  and  the  average  excess  free  energy  per  particle  which  is  simply  replaced 
by  the  excess  free  energy  per  particle  of  a  homogeneous  hard-sphere  fluid  of  the  average 
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density  nr(r).  These  two  average  densities  (or  coarse-grained  densities)  are  not  necessarily 
identical. 


i 

We  choose  here  to  use  the  Carnahan-Starling  equation  of  state  in  our  calculations. 
The  corresponding  excess  hard-sphere  free  energy  per  particle  is 


J ro(n)  =  kT 


y(4-3  y) 

(i  -  y )2  1 


(2.13) 


where  y  =  rend3 /6.  The  expression  for  ^Fo(n)  retains  its  simple  form  as  in  Eq.  (2.13)  when 
extended  to  mixtures  of  hard  spheres  of  equal  size  if  we  identify  n  with  the  total  number 
density  of  particles  of  all  species,  i.e.,  n  =  na. 


In  the  expression  for  the  excess  hard-sphere  free  energy  density  functional,  Eq.(2.l2), 
the  coarse-grained  densities,  n“  (r)  and  n£(r),  are  the  weighted-averages  of  the  local  density 
over  certain  spatial  domains  whose  sizes  are  comparable  to  the  particle  volumes.  This  is 

n"(r)  =  |dVna(r>0(r-r';{n}),  (2.14) 

nra(r)  =  J  d3r'na(r')ra(r-r';{a}),  (2.15) 

where  the  weighting  functions,  v  and  r,  are  required  to  satisfy  the  normalization  conditions 

J d3r'  ua( r-r';{n})  =  J d3r' ra( r-r';{n})  =  l,  (2.16) 

since  these  coarse-grained  densities  both  reduce  to  the  bulk  density  in  the  homogeneous 
fluid  limit. 


In  order  to  complete  the  construction  of  the  excess  free  energy  functional,  prescrip¬ 
tions  for  the  weighting  functions  must  be  specified.  Different  prescriptions  define  different 
models  in  free  energy  density  functional  theory.  The  simplest  criterion  used  to  determine 
the  “optimal”  weighting  functions,  u  and  r,  is  that  they  should  reproduce  the  known  re¬ 
sults  of  the  hard-sphere  fluids  in  the  uniform  density  limit.  For  example,  we  can  require 
both  the  Omstein-Zemike  direct  correlation  function  and  its  first  derivative  with  respect 
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to  the  density  to  be  reproduced  [11].  The  first  free  energy  density  functional  model  which 
incorporated  coarse-grained  densities  was  the  generalized  van  der  Waals  (GVDW)  model 
of  Nordholm  and  co-workers  [13].  The  weighting  functions  of  this  modei  are 

i 

va{r  -  r';  {n})  =  <5(r  -  r'),  (2.17) 


and 


ra(r-r';{n}) 


3 

4  rrd3  ’ 


=  0, 


jr  -  r'|  <  d, 
|r  -  r'|  >  d, 


(2.18) 


where  6( r)  is  the  Dirac  delta  function.  The  GVDW  theory  was  found  to  predict  physically 
unrealistic  behavior  of  inhomogeneous  hard  sphere  fluids  at  high  densities  [14.17],  and 
unsatisfactory  results  in  the  case  of  the  RPM  of  the  electrical  double  layer  [18]. 


There  is  another  free  energy  density  functional  model,  the  generalized  hard-rod  model 
(GHRM),  proposed  by  Robledo  and  Varea  [15],  and  by  Fischer  and  Heinbuch  [16].  The 
GHRM  is  the  three-dimensional  generalization  of  the  one- dimensional  hard-rod  model  and 
is  defined  by 


and 


z/a(r-r';{n})  = 


4  n{d/2)‘ 


:6(d/2  -  |r  -  r'|), 


T«(r-r';{n}) 


3 

4ir(d/2 )3  ’ 


=  0, 


|r  —  r'|  <  d/2, 
|r  —  r'|  >  d/2, 


(2.19) 


(2.20) 


The  counterparts  of  Eqs.(2.19)  and  (2.20)  axe  exact  in  the  one-dimensional  case.  It  has  been 
shown  that  the  GHRM  predicted  qualitatively  correct  but  quantitatively  poor  results  for 
hard-sphere  fluids  confined  between  two  hard  walls  [17].  Yet  the  results  of  the  GHRM  for 
the  RPM  of  the  electrical  double  layer  were  found  to  be  satisfactory,  and  even  quantitatively 
accurate  in  many  cases  [19]. 


Improvements  on  the  GVDW  theory  have  been  made  by  Tarazona  [14]  through  in¬ 
troducing  a  density-dependent  weighting  function  r(r  —  r';  {n}).  The  other  weighting 
function  i/( r  —  r7;  {n}),  however,  has  remained  a  Dirac  delta  function.  In  the  TRZ  model. 
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r  is  assumed  to  expand  in  a  power  series  in  the  coarse- grained  density  nr,  and  the  series 
is  truncated  at  the  third  term: 

r(r  -  r';  {n})  =j=  u;(0)(|r  -  r'|)  4-  u>(1)(|r  -  r'|)nr(r)  +  ic(2)(|r  -  r'|)nr(r)2,  (2.21) 

where  a,(,)(r)  (i  =  0,1,2)  are  density-independent  coefficients.  The  coefficients  of  the 
zeroth-  and  first-order  terms,  u/0)(r)  and  u>(1)(r),  were  determined  by  comparison  with 
the  virial  expansion  of  the  direct  correlation  function  of  a  uniform  hard  sphere  fluid.  The 
coefficient  of  the  second-order  term,  tc(2)(r),  was  calculated  through  an  empirical  fit  to 
the  Percus-Yevick  approximation  of  the  direct  correlation  function  [14].  In  addition,  the 
normalization  condition  for  r(r  -  r';  {n}),  Eq.  (2.16),  implies  that  w(0)(r)  is  normalized, 
and  that  u>^(r)  and  u/2*(r)  have  zero  integrals: 

d3r  u/‘^(|r|)  =  1,  *  =  0, 

=  0,  i  =  l,2. 

The  coefficients  u>(,)(r)  obtained  in  this  manner  axe  optimized.  The  results  are 

3 


/ 


(2.22) 


w 


<0>(r)  = 


w 


(1) 


4flrd3  ’ 

=  0, 

2 


r  <  d, 
r  >  d, 


(2.23a) 


(r)  =  0.475-  0.648  +  0.113  Q  , 

=  0.288  (j'j  -  0.924  +  0.764  -  0.187 

=  0, 

5ird3 


r  <  d, 
d  <r  <  2d, 
r  >  2d, 


(2.23  b) 


u/2)(r)  = 


144 


r  \  2 

[6  -  12  (; 

i)  +  5(; 

5) 

=  0, 


r  <  d, 
r  >  d. 


(2.23c) 


One  realizes  immediately  that  the  GVDW  model  is  just  the  zeroth  order  approximation 
to  the  TRZ  model  by  comparing  Eqs.(2.18)  and  (2.23).  The  TRZ  model  combined  with 
the  Camahan-Starling  equation  of  state  was  found  to  predict  the  best  results  among  those 
three  models  we  discussed  for  hard-sphere  fluids  confined  between  two  hard  walls  [17]. 


Extending  the  TRZ  model  to  mixtures  of  hard  spheres  with  arbitary  size  ratio  involves 
large  complexities  [23];  the  power  series  expansion  of  the  weighting  function  r  is  thus 
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truncated  at  the  linear  term  in  density  as  has  been  done  in  Ref.  23.  In  the  case  of  mixtures 

of  hard  spheres  with  the  same  size,  however,  the  TRZ  model  is  much  easier  to  be  extended. 

Because  only  the  particle  diameter  determines  the  finite  size  effects,  regardless  of  the  other 
* 

parameters  such  as  charge,  it  is  reasonable  to  assume  that  for  mixtures  of  hard  spheres  of 
equal  size  the  weighting  functions,  i/a(r  -  r';  {n})  and  rQ(r  -  r';  {n}),  are  independent  of 
the  particle  species  a  and  take  the  forms 

uQ(r  -  r';  {n})  =  6{r  -  r'), 

and 

ra(r  -  r';  {n})  =  u>(0)(|r  -  r'|)  +  u>(1)(|r  -  r'|)]T  nj(r) 

9 

+  ^(2,(|r  -  r'|)  flS(pWr), 

9y 

where  the  coefficients  are  given  by  Eq.  (2.23).  Eq.  (2.25)  is  essentially  the  same 

as  Eq.  (2.21)  except  that  the  coarse-grained  density  ftr(r)  is  now  identified  with  the  total 
coarse-grained  density  of  particles  of  all  species,  ^(r).  It  follows  from  Eqs.  (2.15)  and 
(2.25)  that 

«;(r)  =  ftL0)(r)  +  <1)(r)nr(  r)  +  ft<,2)(r)ftr(r)2,  (2.26) 

where 

«L°(r)  =  J  <Pr' nQ{r')w(x){\r  -  r'|),  *  =  0,1,2,  (2.27) 

and 

”r(r)  =  n5(r)-  (2.28) 

9 

Summing  over  all  the  equations  of  (2.26)  for  all  species  a  yields 


(2.24) 


(2.25) 


nr(r)  =  ft(0)(r)  +  fta)(r)fT(r)  +  ft(2)(r)ftr(r)2,  (2.29) 


where 

ft«‘'(r)  =  Y,  *?(*)■  (2.30) 

9 
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Eq.  (2.29)  can  then  be  solved  for  nr(r)  to  give 


nr(r)  =  2fiw(r)  j  (l  -  n{1)( r))  +  (l  -  n(1)(r))2  -  4fi(0)(r)n(2)(r) 


-l 


(2.31) 


Eqs.  (2.26)  and  (2.31)  together  determine  the  coarse-grained  density  n£( r),  and  it  is  clear 
that  they  reduce  to  Tarazona’s  results  [13]  in  the  case  of  single-component  hard  sphere 
system. 

It  is  worthwhile  deriving  the  functional  derivatives  of  the  coarse-grained  densities 
with  respect  to  the  local  density  for  the  TRZ  model.  The  results  are 

6n£(r') 


<5na(r) 


=  Sa^6(r'  -  r), 


(2.32) 


6nr(r') 


r(r'  —  r;  nr(r')) 


and 


<§nQ(r)  1  —  n^^r')  —  2n(2)(r')nr(r')  ’ 

6nTg(r')  6flr(  r') 


(2.33) 


-  2ftr(r')  -  «SV) 


(2.34) 


where  SQg  is  the  Kronecker  delta  function. 

Combining  together  all  the  contributions  to  the  Helmholtz  free  energy  density  func¬ 
tional  F({n})  in  Eq.(2.9)  yields 

F({n})  =  J  d3rno(r)  (A3na(r))  -  l] 

a  J 

+  Y1  j  d3rna(r)va(r  ) 


+  T  f  d3rK{ r)F0(fir(r)) 


^  J 


(2.35) 


+  iE/  f  <Prd3r,nUr)nS(*')-^Z!p-[ 

~\kT  J  J  d3rd3r'na(r)ng(r')Acag(\r  -  r'|). 
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The  chemical  potential  of  ions  of  species  a  in  the  electrical  double  layer  is  computed  from 
Eq.  (2.35)  as  the  functional  derivative  of  F  with  respect  to  the  density  of  species  a.  This 
gives 


Ma  =  kT  In  (A3  nQ(r))  +  qQv{r)  4-  u*J(r) 

,3  ,  6n3(r') 


6na{r) 


^o(nr(r')) 


+  £/ 

3  J 

+?^Vn5(r,)["r 6nM  K(r')  j 

-kT^  j  dVn*(r')Aca*(|r-r'|), 


<5n;(r')  ^0(flr(r')) 


where  V(r)  is  the  mean  electrostatic  potential  defined  by 


(2.36) 


0(r)  =  uc(x)+  [  d3r'  ^2 

J  a 


gana(r) 

e|r-r't' 


(2.37) 


Poisson’s  equation  can  be  solved  with  the  appropriate  boundary  conditions  for  the  planar 
electrical  double  layer  to  give  the  mean  electrostatic  potential  in  terms  of  the  ion  density 
distribution  functions  [5].  The  result  is 


0(x) 


)^2gana(x'). 


(2.38) 


Evaluating  the  chemical  potential  in  the  coexisting  homogeneous  bulk  electrolyte  and 
combining  this  with  Eq.  (2.36),  we  obtain  an  integral  equation  for  n0(r)  as  a  function  of 
temperature  and  bulk  ion  densities  n°a .  Because  of  the  planar  symmetry  of  the  problem,  the 
ion  densities  vary  only  in  the  direction  normal  to  the  surface.  Integrating  i/a(r  —  r';  {n}), 
ra(r  -  r';  {n}),  and  Ac0fl(|r  -  i^l)  over  the  directions  parallel  to  the  wall,  i.e.,  the  y-z 
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plane,  leads  to  the  one-dimensional  integral  equation  for  n0(x): 


na(x)  j  qaMx) 

— —  =  exp{  -■ 


kT 


6nuJx')  ^n(n°) 

ona(x)  kT 

6nry(x')  6f0(hr(x')) 


J  dx'na(x') 


n  dT§  (n°) 


Sna(x)  6nry(x')  J 
+  J  dx'  (ng(x')  -  ng)  Acaa(i  -  x')|. 


kT  &nl 

for  x  >  d/ 2,  and 

n0(x)  =  0, 

for  x  <  d/2,  where  n°  is  the  ion  number  density  in  the  bulk  phase  and  n° 


:2.39a) 


(2.396) 

Han°3- 


III.  RESULTS 


The  restricted  primitive  model  of  the  electrical  double  layer,  Eq.  (2.39),  was  solved 
by  means  of  finite  element  techniques  described  elsewhere  [24]  to  obtain  the  equilibrium 
ion  density  distributions  of  1:1  and  2:2  electrolytes.  The  mean  electrostatic  potential 
function  then  followed  immediately  from  Eq.  (2.38).  In  our  method,  the  domain  of  interest, 
d/2  <  x  <  L  {L  is  the  cutoff  value  for  all  the  integrals  in  Eq.  (2.39)),  was  discretized 
uniformly  and  the  trapezoidal  rule  was  used  to  perform  all  the  integrations.  Eq.  (2.39) 
then  becomes  a  system  of  nonlinear  algebraic  equations  for  the  nodal  values  of  the  reduced 
density  gat  =  na(xj)/n° .  We  used  Newton’s  iterative  method  to  solve  for  the  unknowns. 
The  iterations  continued  until  the  Euclidean  norm  of  the  updates  was  less  than  10-10: 

2 N 

where  N  is  the  number  of  nodes  in  the  domain.  Three  to  seven  iterations  were  usually 
required  to  converge  to  a  solution. 


3 


<  10 


-10 


(3.1) 
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A  family  of  solutions  was  composed  of  several  steps  in  the  surface  charge  density  a  at 

a  constant  concentration  c.  For  the  lowest  a,  the  MGC  results  was  used  as  the  initial  guess 

to  the  solution,  whereas  for  the  higher  <r,  the  initialization  was  chosen  by  the  method  of 
» 

parametric  continuation  [25].  The  domain  length  we  used  in  the  calculations  depended  on 
the  concentration  and  the  valence  type  of  the  bulk  electrolyte;  approximately  10  times  of 
the  Debye  screening  length  /c-1  («2  =  ($ir/ekT)  g2  )  was  usually  used.  The  criterion 

employed  to  determine  the  upper  limit  of  the  domain,  £,  is  that  the  ion  densities  at  x  >  L 
should  not  differ  from  their  bulk  values  by  more  than  0.01%.  The  number  of  nodes  N  also 
varied  from  one  concentration  to  another.  A  typical  value  of  N  was  250,  which  was  found 
adequate  in  the  sense  that  further  refinement  did  not  change  the  solution  by  more  than 
1%.  With  241  nodes,  for  example,  our  program  required  less  than  40s  of  Cray-2  CPU  time 
to  obtain  a  solution. 

The  accuracy  of  the  numerical  method  can  be  measured  by  the  degree  to  which  the 
overall  electroneutrality  condition  of  the  interfacial  region, 

roo 

/  dx'  Y]qana(x')  +  <r  =  0,  (3.2) 

Jo  ~ 

is  satisfied.  Eq.  (3.2)  was  satisfied  to  at  least  five  significant  figures,  except  at  the  lowest 
concentrations  considered  (0.01M  for  1:1  electrolytes  and  0.05M  for  2:2  electrolytes),  where 
it  was  satisfied  only  to  four  significant  figures. 

We  used  the  same  conditions  of  calculations  as  those  of  the  MC  simulations  [4]:  d  is 
4.25A,  e  is  78.5,  and  T  is  298K.  The  plasma  parameter  is  T*  =  e2 /ekTd  —  1.6809  in  this 
case.  Dimensionless  quantities  were  used  for  convenience  as  follows.  The  length  is  in  units 
of  d ,  the  surface  charge  density  in  units  of  e/d2,  a*  =  trd^/e,  and  the  mean  electrostatic 
potential  in  units  of  kT/e,  0*(x)  =  t^l(x)/kT,  where  e  is  the  magnitude  of  the  electronic 
charge. 

One  of  the  most  important  physical  quantities  in  the  electrical  double  layer  is  the 
diffuse  layer  potential,  0(0),  which  is  the  potential  difference  between  the  closest  approach 
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(3.3) 


to  the  charged  surface,  d/2,  and  infinity.  It  follows  from  Eq.  (2.38)  that 

0(0)  =  -—/  dx' x'^2qQna(x'). 

6  do  Q 

* 

The  position  of  the  charged  surface  is  shifted  to  x  =  —d/2  in  Eq.  (3.3)  and  in  all  the  results 
to  be  presented  hereinafter.  The  diffuse  layer  potential  can  be  interpreted  as  a  measure  of 
the  thickness  of  the  electrical  double  layer  region  in  which  the  total  charge  qanQ(x)  is 
locally  non-zero. 

1.  1:1  electrolytes 

We  made  calculations  for  1:1  electrolytes  at  concentrations  c  ranging  from  0.01  to  2M 
and  surface  charge  densities  <r*  from  0.05  to  0.9.  In  Table  II  we  report  the  dimensionless 
diffuse  layer  potential  0*(O)  calculated  from  the  GVDW  and  the  TRZ  models,  and  present 
the  MC  data  [4]  and  the  results  from  the  GHRM  [19]  for  comparison.  We  also  list  the  results 
of  the  MGC,  MPB5  [7],  and  YBG  [8]  theories,  and  the  OZ/LMBW  results  by  Plischke  and 
Henderson  [10].  At  low  concentrations  and  low  surface  charge  densities,  the  MPB5  theory 
gives  values  of  0*( 0)  which  agree  best  with  the  MC  results.  The  MPB5  results  for  high  c 
or  a*  are  not  available.  On  the  other  hand,  the  YBG,  the  OZ/LMBW,  and  the  density 
functional  theory  all  give  results  of  comparable  accuracies  at  the  concentrations  and  charge 
densities  we  investigated. 

In  order  to  compare  the  three  different  models  of  the  free  energy  density  functional 
theory,  the  values  of  i/>*(0)  of  these  models  as  a  function  of  <r*  are  plotted  in  Fig.  1.  At 
low  concentrations  and  surface  charge  densities,  our  results  differ  from  the  MC  simulations 
by  only  a  few  percent,  which  is  about  the  level  of  statistical  uncertainty  in  the  MC  data. 
At  high  concentrations  and  surface  charge  densities,  however,  our  values  of  i/>*(0)  are  con¬ 
sistently  lower  them  the  MC  results.  In  this  case  the  TRZ  model  is  obviously  superior  to 
the  GVDW  model  but  poorer  than  the  GHRM;  the  differences  between  the  results  and 
the  MC  values  are  1-5%,  10-20%,  and  20-30%  for  the  GHRM,  the  TRZ  model,  and  the 
GVDW  model,  respectively.  This  is  contrary  to  the  finding  [17]  that  the  TRZ  model  is 
in  general  a  better  approximation  for  inhomogeneous  hard-sphere  fluids  than  the  GHRM. 
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and  confirms  the  suggestion  [19]  that  there  is  some  kind  of  fortunate  cancellation  of  er¬ 
rors  when  the  GHRM  is  combined  with  the  MSA  solution  to  the  electrostatic  correlation 
function  of  the  homogeneous  bulk  electrolytes.  We  should  expect  better  results  from  the 
TRZ  model  than  from  the  GHRM  when  an  improved  treatment  of  the  ion-ion  electrostatic 
interaction  is  used,  since  the  main  approximation  in  our  calculations  is  the  replacement  of 
the  electrostatic  part  of  the  direct  correlation  function  of  the  locally  non-neutral  inhomo¬ 
geneous  electrolyte  by  that  of  the  homogeneous  neutral  bulk  electrolyte.  When  a  neutral 
homogeneous  electrolyte  is  brought  next  to  a  positively  charged  surface,  for  example,  there 
is  a  surplus  of  the  anions  in  the  interface.  The  electrolyte  in  this  region  is  thus  non-neutral 
and  has  different  electrostatic  correlations  from  those  of  a  neutral  electrolyte.  It  may  be 
possible  to  improve  our  results  by  using  the  MSA  solution  to  the  electrostatic  correlation 
function  of  the  homogeneous  but  locally  non-neutral  electrolytes  which  have  the  ion  com¬ 
positions  corresponding  respectively  to  those  found  locally  in  the  double  layer  region.  This 
kind  of  calculation  was  reported  by  Forstmann  and  coworkers  in  the  so-called  local  MSA 
(MSAL)  approach  [9]. 

In  Figs.  2  and  3  we  present  the  reduced  ion  density  na( i)/n°  and  the  dimension¬ 
less  mean  electrostatic  potential  0*(x)  =  c^(x)/kT  profiles  for  c=0.lM  and  cr*=0.3.  At 
this  relatively  low  concentration  and  surface  charge  density,  the  profiles  predicted  by  the 
GVDW  and  the  TRZ  models  are  almost  indistinguishable  (the  difference  is  smaller  than 
5%)  and  in  good  agreement  with  the  MC  predictions  (10%  high).  The  density  and  the 
potential  profiles  are  monotonic  in  this  case. 

One  of  the  most  distinctive  features  of  the  density  profiles  of  the  electrical  double 
layer  is  the  layering  phenomenon  of  the  counterions  near  the  electrode  plate  at  sufficiently 
high  surface  charges.  The  layering  effect  is  caused  by  the  hard-core  exclusion  interplaying 
with  the  purely  electrostatic  response  of  the  ions  to  the  electric  field.  We  show  in  Figs.  4, 
5,  and  6  the  ion  density  profiles  for  1:1  electrolytes  at  c=lM  and  <7*  =0.42,  0.55,  and  0.7, 
respectively.  The  development  of  the  second  layer  of  the  counterions  near  x  =  d  is  plain 
as  ( 7 *  increases.  We  see  in  Fig.  6  that  both  the  TRZ  and  the  GVDW  models  predict 
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approximately  the  correct  position  of  the  second  peak.  In  the  GHRM,  however,  the  peak 
is  shifted  closer  to  the  wall.  The  magnitude  of  the  second  layer  is  not  predicted  well  by  any 
of  the  three  models;  it  is  severely  underestimated  by  the  GVDW  model  and  overestimated 

i 

by  the  GHRM. 

The  mean  electrostatic  potential  profile  in  the  same  conditions  as  in  Fig.  6  is  pie” red 
in  Fig.  7.  In  this  case,  the  MC  simulations  seem  to  predict  a  minimum  in  the  potential, 
but  the  mininum  is  too  shallow  to  be  resolved  from  the  statistical  noise  in  the  data.  Our 
profiles  are  still  monotonic  and  the  deviation  from  the  MC  results  is  the  largest  at  x  =  0. 
that  is,  15%  for  the  TRZ  model  and  35%  for  the  GVDW  model. 

At  even  higher  concentrations  (c=2M,  for  example),  the  MC  simulations  begin  to 
show  the  interesting  effect  of  charge  inversion,  i.e.,  the  coion  density  obviously  exceeds 
the  counterion  density  within  a  certain  region  of  the  interface.  The  result  is  essentially  a 
dipole  layer.  In  Fig.  8  we  present  the  ion  density  profiles  for  c=2M  and  <7* =0.396.  Charge 
inversion  is  seen  in  our  profiles  between  x  =  1.5 d  and  3d  approximately.  In  this  case,  the 
TRZ  model  predicts  results  which  differ  from  the  MC  simulations  by  no  more  than  10%. 
As  a  direct  consequence  of  the  charge  inversion  effect,  the  corresponding  mean  electrostatic 
potential  profile  shown  in  Fig.  9  presentes  an  oscillation  in  the  potential.  The  minimum 
predicted  by  the  TRZ  model  is  —0.166,  which  is  approximately  0.14  more  negative  than 
that  in  the  MC  profile.  Compared  with  the  MC  results,  the  position  of  our  minimum  is 
shifted  toward  the  charged  wall. 

2.  2:2  electrolytes 

The  2:2  electrolyte  is  known  to  be  a  more  stringent  test  of  the  electrical  double  layer 
theory  because  the  ionic  interactions  are  much  stronger  than  those  in  the  1:1  electrolyte. 
We  computed  the  density  and  the  mean  electrostatic  potential  profiles  for  c=0.05  and 
0.5M,  and  <7*  varying  from  0.025  to  0.4.  Some  selected  values  of  t/>*(0)  are  also  displayed 
in  Table  II  to  compare  with  the  MC  data  and  results  of  other  theories.  Again,  we  plot 
0*(O)  as  a  function  of  <7*  in  Fig.  10.  Our  results  show  a  distinct  maximum  in  the  diffuse 
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layer  potential,  as  predicted  by  the  MC  simulations;  but  they  agree  only  qualitatively  with 

the  MC  data.  The  disagreement  is  as  high  as  14%  at  Q.05M  and  28%  at  0.5M.  The  TRZ 

and  the  GVDW  models  again  give  indistinguishable  results  at  low  <7*. 

* 

The  density  profiles  of  a  2:2  electrolyte  at  c=0.5M  and  a*  =0.1704  axe  shown  in 
Fig.  11.  There  is  charge  inversion  at  this  moderate  concentration  and  surface  charge 
density.  It  appears  that  charge  inversion  is  easier  to  get  with  2:2  electrolytes  than  with 
1:1  electrolytes  because  of  the  stronger  interactions  between  the  ions  and  between  the  ion 
and  the  charged  surface  in  divalent  electrolytes. 

The  mean  electrostatic  potential  profile  at  the  same  conditions  as  in  Fig.  11  is  pre¬ 
sented  in  Fig.  12.  The  potential  is  oscillatory  with  a  minimum  of  —0.22  at  x  =  0.5d 
approximately.  The  minimum  in  our  profile  is  18%  lower  than  that  in  the  MC  results. 


IV.  SUMMARY 


Free  energy  density  functional  theory  has  been  applied  to  the  restricted  primitive 
model  (RPM)  of  the  electrical  double  layer.  The  hard-core  repulsive  contribution  to  the  free 
energy  was  incorporated  by  a  nonlocal  excess  hard-sphere  free  energy  density  functional 
constructed  on  the  generalized  van  der  Waals  (GVDW)  model  or  the  Tarazona  (TRZ) 
model.  The  electrostatic  part  of  the  ion-ion  correlation  function  of  the  electrolyte  in  the 
interfacial  region  was  approximated  by  the  mean  spherical  approximation  (MSA)  to  that 
of  the  homogeneous  bulk  neutral  electrolytes.  For  each  model,  the  calculations  required  a 
hard-sphere  equation  of  state  as  input.  The  Camahan-Starling  equation  of  state  was  used. 

The  free  energy  density  functional  theory  predicts  both  the  ion  density  and  the 
mean  electrostatic  potential  profiles  that  are  in  good  agreement  with  those  of  the  Monte 
Carlo  (MC)  simulations.  The  theory  is  able  to  predict  correctly  the  layering  effect  of 
the  counterions  near  the  highly  charged  surface  and  the  charge  inversion  phenomenon  at 
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sufficiently  high  concentrations  and  surface  charge  densities.  At  low  concentrations  and 
surface  charge  densities,  the  GVDW  and  the  TRZ  models  give  almost  identical  results. 
As  either  the  concentration  or  the  surface  charge  density  increases,  the  TRZ  model  is 

i 

obviously  superior  to  the  GVDW  model.  The  diffuse  layer  potentials  predicted  by  the 
TRZ  model  are  poorer  than  those  by  the  generalized  hard-rod  model  (CHRM)  reported 
in  Ref.  19.  We  expect  that  the  theory  will  be  improved  by  using  the  MSA  solution  to  the 
electrostatic  correlation  function  of  the  homogeneous  non-neutral  electrolytes  which  have 
the  same  ion  compositions,  respectively,  as  those  found  locally  in  the  interface.  This  is  a 
matter  for  future  investigations. 

The  relative  accuracy  of  the  different  theories  of  the  RPM  of  the  electrical  double 
layer  is  mainly  of  theoretical  interest  [6],  since  the  RPM  is  such  a  greatly  simplified  model. 
Of  far  greater  relevance  is  to  extend  these  different  theories  to  more  realistic  models  that 
take  into  account  more  and  more  physical  features  of  the  electrical  double  layer.  These 
features  include  the  solvent  effects  (electrostriction,  dielectric  saturation),  the  ion  effects 
(ion  polarizability,  ion-concentration  dependence  of  the  solvent  dielectric  constant),  the 
asymmetric  electrolytes  (size  asymmetry,  charge  asymmetry),  and  the  surface  properties 
(discrete  charge  distributions,  image  forces,  non-planar  surfaces),  etc.  The  simplest  model 
which  incorporates  the  solvent  effects  is  the  hard  sphere  ion-dipole  mixtures  (HSIDM) 
model  [5,6],  in  which  the  solvent  is  represented  as  a  fluid  of  hard  spheres  with  embedded 
point  dipoles.  Although  this  is  still  a  crude  model  of  an  electrolyte,  it  is  certainly  an 
improvement  over  the  primitive  model.  The  extension  of  the  density  functional  theory  to 
the  HSIDM  model  is  straightforward.  It  requires  an  extra  term  to  incorporate  the  finite 
size  effects  of  the  solvent  molecules  and  two  extra  terms  to  include  the  ion-dipole  and 
dipole-dipole  interactions  in  the  expression  of  the  free  energy  density  functional.  We  are 
currently  working  on  this  model. 
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FIGURE  CAPTIONS 


Fig.  1  Dimensionless  diffuse  layer  potential,  ev{0)/kT,  as  a  function  of  the  surface  charge 
density,  ad2 /e,  for  1:1  electrolytes  at  c=0.01,  0.1,  1,  and  2M.  The  solid  curves,  the 
dot  curves  and  the  dashed  curves  are  the  results  of  the  TRZ,  the  GVDW  and  the 
GHRM  models,  respectively.  The  dots  and  the  squares  represent  the  MC  data. 

Fig.  2  Reduced  ion  density  profiles,  n(x)/n°,  for  a  1:1  electrolyte  at  c  =  0.1M  and  a’  =  0.3. 
The  solid  curves,  the  dot  curves  and  the  dashed  curves  are  the  results  of  the  TRZ. 
the  GVDW  and  the  GHRM  models,  respectively.  The  dots  represent  the  MC  results. 
The  charged  wall  is  at  x  =  —0.5 d. 

Fig.  3  Dimensionless  mean  electrostatic  potential  profile,  eip(x)/kT,  for  a  1:1  electrolyte  at 
c  =  0.1M  and  <7*  =  0.3.  All  symbols  as  in  Fig.  2. 

Fig.  4  Reduced  ion  density  profiles,  n(x)/n°,  for  a  1:1  electrolyte  at  c  =  1M  and  a “  =  0.42. 
All  symbols  as  in  Fig.  2. 

Fig.  5  Reduced  ion  density  profiles,  n(x)/n°,  for  a  1:1  electrolyte  at  c  =  lM  and  cr*  =  0.55. 
All  symbols  as  in  Fig.  2. 

Fig.  6  Reduced  ion  density  profiles,  n(x)/n°,  for  a  1:1  electrolyte  at  c  =  lM  and  a*  —  0.70. 
All  symbols  as  in  Fig.  2. 

Fig.  7  Dimensionless  mean  electrostatic  potential  profile,  exjj(x)/kT,  for  a  1:1  electrolyte  at 
c  =  lM  and  a *  =  0.70.  All  symbols  as  in  Fig.  2. 

Fig.  8  Reduced  ion  density  profiles,  n(x)/n°,  for  a  1:1  electrolyte  at  c  =  2M  and  a*  —  0.396. 
All  symbols  as  in  Fig.  2. 

Fig.  9  Dimensionless  mean  electrostatic  potential  profile,  ei l>(x)/kT,  for  a  1:1  electrolyte  at 
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c  =  2M  and  a’  =  0.396.  All  symbols  as  in  Fig.  2. 


Fig.  10  Dimensionless  diffuse  layer  potential,  erp(0)/kT ,  as  a  function  of  the  surface  charge 
density,  acP /e, 'for  2:2  electrolytes  at  c=0.05  and  0.5M.  All  symbols  as  in  Fig.  1. 

Fig.  11  Reduced  ion  density  profiles,  n(x)/n°,  for  a  2:2  electrolyte  at  c  =  0.5M  and  a*  = 
0.1704.  All  symbols  as  in  Fig.  2. 

Fig.  12  Dimensionless  mean  electrostatic  potential  profile,  txp(x)/kT ,  for  a  2:2  electrolyte  at 
c  =  0.5M  and  cr*  =  0.1704.  All  symbols  as  in  Fig.  2. 
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Table  I.  Theories  of  the  RPM  of  the  electrical  double  laver 


Integral  Equation  Methods 


Density  Functional  Theories 


Table  II.  Diffuse  layer  potential  'P'(O)  =  e'I'(0)/fcT 


c 

<7* 

MC“ 

MGC 

BGY* 

’  MPB5C 

PH* 

GHRMe 

GVDW  f 

TRZ' 

1:1  electrolytes 

0.01M 

0.10 

5.05(0.05) 

5.44 

— 

5.08 

4.58 

5.26 

5.242 

5.243 

0.1M 

0.30 

4.63(0.03) 

5.34 

5.00 

4.74 

4.37 

4.76 

4.49 

4.54 

1.0M 

0.10 

1.09(0.06) 

1.40 

1.055 

1.03 

1.06 

1.03 

1.00 

1.01 

0.25 

2.13(0.05) 

2.79 

2.31 

2.10 

2.22 

2.18 

2.01 

2.05 

0.42 

3.08(0.10) 

3.74 

3.46 

3.02 

3.23 

3.23 

2.68 

2.55 

0.55 

4.15(0.15) 

4.26 

4.21 

— 

4.22 

4.12 

3.12 

3.61 

0.60 

4.38(0.11) 

4.43 

4.48 

— 

4.68 

4.52 

3.29 

3.96 

0.70 

5.71(0.14) 

4.74 

5.02 

— 

5.76 

5.41 

3.69 

4.79 

2.0M 

0.396 

2.29(0.09) 

2.99 

2.303 

— 

2.29 

2.19 

1.68 

1.55 

2:2  electrolytes 

0.05M 

0.20 

1.33(0.02) 

2.61 

1.81 

1.36 

1.18 

1.59 

1.539 

1.543 

0.5M 

0.1704 

0.63(0.04) 

1.36 

0.638 

0.537 

0.69 

0.57 

0.528 

0.532 

».  G.  M.  Torrie  and  J.  P.  Valleau,  Ref.  4.  Statistical  uncertainty  is  shown  in  parenthesis. 

b.  C.  Caccamo,  G.  Pizzimenti  and  L.  Blum,  Ref.  8;  1986. 

c.  C.  W.  Outhwaite  and  L.  B.  Bhuiyan,  Ref.  7;  1986. 

d.  M.  Plischke  and  D.  Henderson,  Ref.  10;  1988. 

e.  L.  Mier-y-Teran,  S.  H.  Suh,  H.  S.  White,  and  H.  T.  Davis,  Ref.  19. 

f.  This  work.  Results  of  the  GVDW  model. 

g.  This  work.  Results  of  the  TRZ  model. 


2  4  6  0  10 

x/d 


Figure  2 


(x)/n 


Figure  4 


Figure  6 


Figure  7 


3 


2.5 


Figure  9 


/3ef(0) 


Figure  10 


Figure  11 


